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O ■ ABSTRACT 

C^ I By analysing models of the young massive cluster R136 in 30 Doradus, set-up using 

O !■ the herewith introduced and publicly made available code McLuster, we investigate 

and compare different methods for detecting and quantifying mass segregation and 
substructure in non-seeing limited iV-body data. For this purpose we generate star 
cluster models with different degrees of mass segregation and fractal substructure and 
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~* ' analyse them. 
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We quantify mass segregation by measuring, from the projected 2d model data, 
the mass function slope in radial annuli, by looking for colour gradients in radial 
\^ colour profiles, by measuring Allison's A parameter, and by determining the local stel- 

lar surface density around each star. We find that these methods for quantifying mass 
segregation often produce ambiguous results. Most reliable for detecting mass segre- 
gation is the mass function slope method, whereas the colour gradient method is the 
pS 1 least practical in an R136-like configuration. The other two methods are more sensitive 

to low degrees of mass segregation but are computationally much more demanding. 
We also discuss the effect of binaries on these measures. 

Moreover, we quantify substructure by looking at the projected radial stellar den- 
sity profile, by comparing projected azimuthal stellar density profiles, and by deter- 
mining Cartwright & Whitworth's Q parameter. We find that only high degrees of sub- 
structure affect the projected radial density profile, whereas the projected azimuthal 
density profile is very sensitive to substructure. The Q parameter is also sensitive to 
substructure but its absolute value shows a dependence on the radial density gradient 
C$ ' of the cluster and is strongly influenced by binaries. 

Thus, in terms of applicability and comparability for large sets of -ZV-body data, 
the mass function slope method and the azimuthal density profile method seem to 
be the best choices for quantifying the degree of mass segregation and substructure, 
respectively. The other methods are computationally too demanding to be practically 
feasible for large data sets. 

Key words: galaxies: star clusters: individual: R136 — galaxies: star formation — 
methods: data analysis — Magellanic Clouds 



1 INTRODUCTION cepted picture, star clusters form in three stages: first, a 

cold molecular cloud collapses and forms stars along fila- 
Understandme the process or star cluster formation is vital ,, , ,,. „ „ ,, ,, „ 

, r ments throughout this collapse. Secondly, the massive O- 

for astrophysics since most, if not all stars are born in a , ,-, ,. ,. „ ,, • , , , ., , 

17 — : — - — 7 — '7-11 -1. and B-stars start radiating oft the residual gas until only 

clustered mode (Lada & Lada 2003). In the commonly ac- , , , , , „ , r . „ , 

' ' a more or less bound ensemble of stars is left, subse- 

quently, the newly formed star cluster evolves dynamically 
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the birth process and the duration of the subsequent dy- 
namical dissolution depends crucially on the star formation 
efficiency, i.e. how much of the cold gas gets transformed 
into stars. Moreover, it depends on the structure of the gas 
cloud and the distribution of the forming stars. That is, for 
a more massive ensemble of stars the gas expulsion process 
will be more violent, whereas for low-mass configurations gas 
expu ls ion will happen more ad i abatically (IGever fc Burkert 
200 it iKroupa fc Boilvl 12002); iGoodwin fc Bastianl 12006 



For exampl e lLada fc Ladal d2003h. iTeixeira et al.l J2006h . 
lAllen et all (|2007h as well as iGennaro et al.l (|201lh find 



Basti an fc Goodwin! 20061; Baumgar dt fc Kroupal 12007. 
Baumgardt. Kroupa fc Parmentierl 120081 ). A further influ- 



ence on the survival rate is given by the degree of mass 
segregation, i.e. when the most massive stars of an ensemble 
are located deeper within the forming cluster, they will have 
a more destructive influence on the subsequent evolution 
of th e star cluster (jVesperini. McMillan fc Portegies ZwartJ 
2009). A yet open question in this respect is whether or not 
star clusters form with primordial mass segregation, or if the 
initial mass function (IMF) of stars is the same throughout 
the whole star forming complex. 

Observations of young clusters indeed suggest the cxis- 
tence of primordial mas s segregation fe.g. [Stolte et al. 2002; 
iBontemps et al.l |2010bJ) . But such observed mass segrega- 
tion is not necessarily due to variations of the IMF, as mass 
segregation can also develop quickly during the first few 
100,000 years of cluster formation through dynamical re- 
laxation. The timescale of this process is proportional to 
the relaxation time of th e configuration and the mass ra- 
tio of the forming stars (|Spitzerl 1 19871 ). Moreover, recent 
investigations show that initially substructured configura- 
tions can de velop mass segregation on significantly s horter 
timescales ( McMillan. Vesperini fc Portegies Zwartl 1 2007; 



Allison et al]|2009t iMoeckel fc Bonnelll 120091 ; lAllison et ail 
2010l ; lYu. de Griis fc Cherj|201ll ). In this picture, (fractal) 



substructures, which may have much shorter relaxation 
times, can segregate before they merge to form the final 
star cluster. Hence, substructure in young star clusters is 
not only a sign of a system not being virialised, but may 
also play a vital role in the process of star cluster formation. 
Thus, two major aspects of young star clusters have to be 
investigated in detail at the different stages of star cluster 
formation: the degree of mass segregation and the degree of 
substructure. 

These questions may be addressed by means of col- 
lisional iV-body computations. But such numerical inves- 
tigations also require a choice of initial conditions. The 
question therefore remains what configuration star clus- 
ters have at the stage where dynamical investigations can 
set in. From hydrodynamical computations of collapsing 
gas clouds it appears that star formation may be hap- 
pening in a hierarchical, f ractal fashion (IKlessenl l200lt 
IBonnell, Bate, fc Vindlioolil ; iBonnell fc Batell2006h . During 
the subsequent dynamical evolution this substructure is 
erased and i ndeed a significant degree of mass segregation is 
established (jMaschberger et alj|2010r ). But hydrodynamical 
computations only reach gas cloud masses of a few 10 3 Mq, 
thus they do not shed any light on star formation in star- 
burst regions or on the formation of globular clusters, nor 
can they account for the self-regulation induced by stellar 
feedback, yet. 

Observations of young embedded star clusters show 
a similar picture like hydrodynamical computations. 



many young clust ers to show substructure and to be asym- 
metric. Moreover. iGutermuth et al . (2005, 2008) find young 
embedded clusters in near-IR data to be azimuthally asym- 
metric with a high degree of substructure. The star forma- 
tion sites appear filamentary and elongated over scales of 
several parsec. Such observations suggest that young clusters 
expand, get more s ymmetric and lose subs t ructure with on- 
going gas removal (IGutermuth et alj|2008t IBontemps et al.l 
l2010aT ). The nearby, more evolved Orion Nebula Cluster, 
for example, shows a high degree of mass segregation which 
appears to be inconsistent with its current re laxation time 
l|Hillenbrand fc Hartma"nnlll998t lKroupall2002l 'l. In this pic- 
ture, this may indicate that the ONC also formed with a 
high degree of substructure. 

But does this picture also apply to starbursts, and to 
the formation of globular cluster-like objects? Massive star 
formation sites with stellar masses above 10 4 Mq in the 
Milky Way like NGC 3603, Westerlund 1 and the Arches 
cluster, are rare and mostly heavily obscured by interstel- 
lar dust. Nevertheless, mass segregation as well as high 
degrees of substructure have been r eported for those ob- 
jects (e.g. [Stolte et all |2002| , I2006J ; iBrandner et~al] 120081 : 
IGennaro et al-lbOHl ). 

The most massive star formation site in the Local 
Group is the 30 Doradus complex in the Large Magellanic 
Cloud. It is the only nearby starburst region, which makes 
it the "Rosetta Stone" for unde rstanding event s such as the 
formation of globular clusters (|Walbornlll99lr ). This paper 
therefore addresses this star formation site and aims at in- 
vestigating mass segregation and substructure in this com- 
plex. For this purpose, we create models of the young mas- 
sive cluster R136 which is forming in the 30 Doradus com- 
plex (Sec. [2] & [3]). We then use various methods from the 
literature (Sec. U} to detect and quantify mass segregation 
and substructure (Sec.[S|. Our aim is to test and calibrate 
these methods in order to be able to apply them to iV-body 
computations which will be presented in a follow-up inves- 
tigation, and to discuss their applicability to observational 
data. Furthermore, Appendix A contains a manual on our 
new star-cluster initialisation code McLuster. 



2 R136 

The young massive cluster R136 is at the heart of the 30 
Doradus (30 Dor) star forming region in the Large Mag- 
ellanic Cloud (LMC). 30 Dor is known as the largest HII 
region in the Local Group with more tha n 8 x 10 5 Mq of 
ionised gas within a radiu s of about 100 pc (|Kennicuttlll984l ; 
iMalumuth fc Head Il99"4r i . The densest region of 30 Dor is 
the central star cluster NGC 2070 with a half-light radius 
of about 22 pc. R136, the centre of this cluster, has always 
posed a challenge to high-resolution observations. This ob- 
ject is crucial for our understanding of star formation, since 
there is no other comparable starburst site in the local Uni- 
verse. But due to its distance, its high central brightness and 
its substantially var ying extinction, R1 36 is hard to access 
observationally fe.g. lBrandl et alj|l996h . 

R136 was found to be about 3 My r old, but its con- 
stituent stars show some age spread IjBosch et al.l l200ll ; 
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lAndersen et al. 2009). Due to the absence of red supergiants 
and due to the presence of several Wolf-Rayet stars, the age 
of the oldest population of stars within R136 can be limited 
to 3-5 Myr (jBrandl et alj|l99q ). Some O stars in the cen- 
tre of R136 may be less than 2 M yr old, though, indicatin g 
a complex star formation history (Ma ssev fc Hunterl ll998). 
The metallicity of the whole 30 Dor co mplex was found to 
be ab out half solar value, i.e. Z ~ 0.01 IjLebouteiller et al.l 
I2008T I. 

Mass estimates for the stellar component of 

R136/NGC2070 ra nge from about 1.4 / 10 4 M Q 

|Malumuth fc Head 
UBV photometry, 



1994), 



to 



estimated from HST 
about 4.5 x 10 5 M 
IjBosch. Terlevich fc Terlevicbl |2009j), estimated from 
multi-epoch stellar velocity dispersion data. Recent high 
precision photometry y ields masses between ~ 5.5 x 10 4 Mp 



llCrowther et all |2010| ). and > 10 s M Q |Andersen et al.1 
2009). 

R136 is an important test bed for our understanding 
of the initial mass function of stars, and thus has been the 
subject of several mass-function investigations. Especially 
the mass function at the high-mass end and a probable upper 
limit of stellar masses have been primary targets of such 
investigations. Down to about 1.1 Mq the mass function of 
R136 was found to be in agreement with a Salp eter slope of 
2.35 l|Hunter et al.1 1 19951 ; lAndersen et"aLll2009n , with some 
of these investigatio ns reporting small radial dependencies 
of the mass function (IBrandl et al.lll996l ; ISelman et al.lll999l ; 
ISirianni et 111120001 ). 

The masses of the brightest objects in R136 are very un- 
certain. Due to resolution limits, R136 was once believed to 
be a single supermassive star of more than 1, 000 Mq. Only 
with Speckle interferometry and the superior resolution of 
the HST it was later found to be a very dense star clus- 
ter (s ee e.g. IWeigelt et al.lll99ir ). Recently. ICrowther et al.l 
(|2010l ) identified four stars in the centre of R136 to have 
masses between 165-320 Mq , again challenging the com- 
monly believed upper ste llar mass limit of about 150 Mq 
l|Weidner fc KroupalEfflol ). 

R136 most proba bly has a high binary frac- 



tion. ISelman et al.l ((1999) and lBosch. Terlevich fc Terlevich! 
(2009) find t hat at least all O- and B-stars are in binaries, 
even though ICrowther et al.l (|2010l ) find that the four most 
massive stars are most likely single stars. But the observed 
velocity dispersion seems to be dominated by binary m otion 
l|Bosch et alj|200ll ; iBosch. Terlevich fc Terlevichll2009h . 

Like many young clusters in the LMC, R136 follows a 
shallow power-law density profile wi thout any visible trun- 
cation at large radii, as specified bv lElson. Fall fc Freeman! 

ma), 



p(R) = p (l + R 2 /a 



2\(— Y/2) 



(1) 



where R is the projected radius, po is the central density, a is 
a scale radius and 7 is the power-law slope. For R136's num- 
ber densit y profile, 7 was fou nd to be about 1.85 for mas- 
sive stars IjHunter et al.lll995T ). Its surface brightness profile 
shows a variety of power-law slopes, depending on the in- 
strument and fi lter which was used. T he 7 values lie between 
1.7 in F336W jCampbell et al.lll992l) and ab out 2 in F555W 
l|Selman et aT]fl999h . I Andersen et al.1 l|2009h found values of 
7 = 1.54 and a = 0.025 pc in F160W, whereas accord- 
ing to iMackev fc Gilmorel (|2003l ) 7 = 2.43 for F555W and 



F814W. [Campbell et ail (|2010h measured 7 = 1.8 in optical 
data (V and I), but 7 = 1.6 i n near-infrared (H and K) im- 
ages. [ Hunter et alj |l995j) and lMcLaughlin fc van der Marell 
f;2005l ) found 7 ~ 2 in F555W. In none of the above investi- 
gations a proper core could be identified, such that the scale 
radius a was in all cases found to be of the order of the 
resolution limit. 

R136 does not appear to be kinematically relaxed 
l|Selman et al.lll999P l. In several studies evidence of dynam- 
ical substructure has been found , such as a ring of mas- 
sive stars at about 2-3 pc radius (|Malumuth fc Heap|ll994J ; 
IBrandl et al. 2007), or a shell of massive stars at a ra- 
dius of 6 pc (|Hunter et al.lll995i ) - structure which should 
be quickly erased by two-body relaxation. Moreover, the 
azimuthal density pro file of R136 shows strong variations 
IjCampbell et al.1 120101 ). Such substructure and asymmetry 
also causes bumps i n the surface brightness profile at the 
corre sponding radii fSel man et al.lll999l ; iMalumuth fc Heap! 
[l994j), making the investigation of R136 difficult. 



3 MODELS 

To calibrate our methods for detecting mass segregation and 
substructure, wegenerated two sets of star cluster models 
with McLusteeQ with properties similar to those of R136 
as seen today (see previous section). In addition, we pro- 
duce one set of mass segregated models without binaries for 
comparison, and two more sets of models with substructure 
using different random seeds for estimating the stochastic 
scatter of the results. 

All calibration models share the same basic properties, 
i.e. the clusters have a total mass of 10 5 M q, consis t ing o f 
about 170,000 stars drawn from the canonical lKroupal (|2001h 
IMF following a power-law index, a, of 1.3 for stellar masses 
between 0.08 Mq and 0.5 Mq, and a — 2.3 for stellar masses 
larger than 0.5 Mq. As the upper stellar mass limit we chose 
100 Mq since stellar evolution models for higher masses are 
rather unreliable and therefore cannot be properly mod- 
elled b y the stellar evolution routines we use. lCrowther et al.l 
((2010j) find 4 stars in R136 to exceed the mass of 165 Mq , 
though. They also estimate the total number of stars with 
initial masses above 100 Mq to be about 14 within a radius 
of 5 pc. While this limitation may only have a negligible 
effect on most measures we test here, it may well have a sig- 
nificant effect on the colour-gradient method to detect mass 
segregation due to their high luminosities (see Sec. |3}. 

The cluster stars have a metallicity of Z = 0.01 and 
they were evolved from the zero-age main sequence to an 
age o f 3 Myr using the SSE routine l|Hurlev, Pols fc Tout! 
120001 ) within McLuster (see Appendix A). Thus, the most 
massive stars have masses of about 80 Mq. This was done in 
order that the cluster stars have comparable colours like the 
stellar population of R136 which is especially important for 
detecting mass segregation with the colour-gradient method. 

The binary fraction, fbin, is 1.0 in all calibration models 
(except, of course, for the set of models without binaries), 



www . astro . uni-bonn . de/~akuepper/mcluster/mc luster . html 
or www. astro .uni-bonn.de/~webaiub/germaii/dowiiloads .php 
A manual on our new publicly available cluster-initialisation code 
McLuster is provided in Appendix A. 



4 A.H.W. Kiipper, Th. Maschberger, P. Kroupa and H. Baumgardt 









■ -ijl- • 


$$- : - 


• 



Figure 1. Logarithmic intensity maps of 6 clusters models which where set-up with McLuSTER. The field-of-view is 20 pc X 20 pc. The 
clusters show different degrees of mass segregation. The mass segregation parameter, S, of the clusters is 0.0 (upper left), 0.4 (upper 
right), 0.7 (middle left), 0.9 (middle right), 0.95 (lower left), and 1.0 (lower right), respectively. All models follow the same mass profile 
and extend out to a radius of 20 pc. 
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Figure 2. Logarithmic intensity maps of 4 cluster models which where set-up with McLuster. The field-of-view is 20 pc X 20 pc. The 
clusters are not mass segregated (5 = 0) and show different degrees of substructure. The fractal dimension, D, of the clusters is 2.6 
(upper left), 2.3 (upper right), 2.0 (lower left), and 1.6 (lower right), respectively. All models follow the same mass profile and extend 
out to a radius of 20 pc. Note that these clusters are not meant to look like real, existing objects but shall rather give a set of models 
with a smoothly increasing degree of substructuredness. 



i.e. all stars are in binaries. The binaries were set up using or- 
dered pairing for stars more massive than 5 M Q following the 
method introduced by Oh & Kroupa (in prep.), i.e. the most 
massive star is in a binary with the second-most massive 
star, the third with the forth, and so on. Stars with masses 
below this threshold were paired randomly. The value of 
5 M© is somewhat arbitrary, but it rests on the observation 
that late-type stars with masses below 1 — 5 Mq follow well 
defined, simple pairing rules (random pairi ng from the IMF , 
the initial period distribution function from lKroupal (|1995bh , 
thermal eccentricity distribution), while massive stars with 
masses larger than about 10 Mq te nd to have sim ilar com- 
ponent masses and shorter periods (|Kroupall201lr ). 

The orbital e l ements of the binaries were generated 
using the iKroupal jl995bl ) pe riod distributi on and a ther- 
mal eccentricity distribution (IKroupal 12008 ). As shown in 
iKiipper. Kroupa fc Baumgardtl IMPOST ), this results in a sig- 



nificant number of binaries which are too wide to be bound 
inside the very dense environment of our models. The mean 
kinetic energy, or 'dynamical temperature', of centre-of-mass 
particles in such a configuration is about 1 km 2 s -2 . Assum- 
ing that all binaries with binding energies lower than this dy- 
namical temperature are unbound or get disrupted quickly 
(Hcggic 1975), the effective binary fraction would be about 
0.5, but about 1.0 among the h igh mass stars which is consis- 
tent with recent o bservations (|Bosch, Terlevich fc Terlevichl 
120091 ) , even though lCrowther et al.l ( 20100 find the four most 
massive stars in R136 not to be in binary systems. But those 
few objects appear to be peculiar in many aspects which may 
well be due to the frequent strong gravitational encounters 
they must experience. Thus, they may be neglected here in 
this respect. 

The density profile, p(R), of the models was cho sen to 
be an EFF profile (eq.fTl lElson. Fall fc Freemanlll987T > with 
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a scale radius, a, of 0.1 pc and a slope of 7 = 2.0. The 
2D density profile was deprojected within McLuster and 
used to generate the 3D cluster configuration (for details see 
Appendix A). The (infinitely extended) EFF profile was cut 
off at a radius of 20 pc, as we are mainly interested in the 
inner ~10 pc. The central density, po, was fixed such that 
the total mass within this radius is 100, 000 Mq. 

To guarantee a good comparability between the differ- 
ent models, we chose the same random seed for all models 
such that the stellar populations in all of them are the same, 
and only the spatial distributions of the stars are different. 
Just the two additional sets of models with substructure 
have each a separate random seed to estimate the effect of 
stochasticity of the initialisation process on the results. 

The two calibration sets of models with binaries differ 
only in one parameter, one has a varying degree of mass 
segregation, and the other has a varying degree of fractal 
substructure. 

(i) We produced 10 models with mass segregation val- 
ues, S, ranging from 0.0 (unsegregated) to 0.9 in steps of 
0.1, and another 10 models with values of S from 0.91 to 
1.0 (completely segregated) in steps of 0.01 (Fig. [TJ. For 
segregating the clusters we chose the me thod suggested by 
iBaumgardt. De Marchi fc Kroupal (|200a ). which preserves 
the desired (mass) density profile. We refer to Appendix A 
for details on how the intermediate steps between unsegre- 
gated and completely segregated clusters were set up with 
McLuster. In short: with a higher value of S, more mas- 
sive stars get higher probabilities to be placed deep in the 
potential well of the cluster, i.e. in the centre. Lower values 
of S correspond to more similar probabilities between high 
and low mass stars, i.e. more random distributions. 

(ii) We generated 3 x 15 models with fractal initial condi- 
tions (Fig. [2}. The fractal dimension, D, of the models was 
varied from 3.0 (non-fractal) to 1.6 (highly fractal) in steps 
of 0.1. Each set of 15 models has a different random number 
seed to measure the stochastic scatter between different re- 
alisations. For fractalizing the stellar distributions, we used 
the procedure described in Appendix A. In short: we set up 
a fractal distribution of sta rs within a unit-sphere follow- 
ing roughly the method of iGoodwin fc Whitworthl (120041 ) 
and "folded" this distribution with the desired density pro- 
file. In this way, we end up with a radially concentrated 
but fractal distribution of stars. This is reminiscent of the 
filamentary and radially oriented structure of dense star- 
forming gas in contracting molecular cloud cores. Moreover, 
this way of producing a radial but fractal distribution is an 
important innovation for testing substructure-detection al- 
gorithms, since young clusters like R136 are neither purely 
radial nor purely fractal (see Sec. [2}. 

Additionally, we produce another set of mass segregated 
models like the one above but without binaries, to determine 
the influence of binary stars on the methods for detecting 
mass segregation. 



4 METHODS 

There have been several attempts to detect and quantify 
mass segregation and substructure in (young massive) star 
clusters. Here we are going to apply some of these techniques 



to our models of R136 to test their feasibility and to get some 
comparability among them. In a follow-up investigation we 
will apply some of these methods to iV-body computations 
to follow the time evolution of mass segregation and sub- 
structure. Since some of the methods are computationally 
intensive when applied to a cluster with 170,000 stars, we 
apply a low mass cut-off at 1.1 Mq for most methods which 
leaves us with about 15,000 stars. This furthermore reflects 
the common situation faced with observational data which 
often suffers from incompleteness and crowding. 



4.1 Mass segregation 

We consider the following methods for the detection and 
quantification of mass segregation: 

(i) In analogy to the work of lAndersen et al.l (|2009j) on 
R136, we measure the mass function slope, a, of stars above 
1.1 Mq at projected radii between 3 pc and 7 pc. In this 
radial range Andersen et al. achieve reasonably high com- 
pleteness levels. For comparison, we also do this for all stars 
in the cluster and for all stars within a projected radius 
of 3 pc from the cluster centre. The slope a we determine 
wi th the Modified M aximum Likelihood estimator (MML) 
of Maschberger & Kroupa (2008). A standard deviation is 
estimated using 100 Monte Carlo subsets of stars. Since all 
models are set up with a mass function slope of 2.3 for stel- 
lar masses above 0.5 Mq, mass segregated clusters should 
show a steeper slope outside 3 pc and a shallower slope in- 
side this radius. The choice of 3 pc is somewhat arbitrary 
but is mean t to establish some com parability to the results 
achieved by lAndersen et al.l |2009). 

(ii) We meas ure radial colour gradie nts following the 
methodology of iGaburov fc Gielesl (|2008T l. From the stellar 
radii and corresponding effec tive temperatures, which are 
provided by the SSE routine dHurlev. Pols fc Tout|[2000h in 
Nbody6 as well as McLuster (see Appendix) , we compute 
B-, V- and /-band m agnitude s for e ach star. We use the 
algorithm described in iFlowerl (1 19961 ) to first compute the 
bolometric correction, BC, and the colour index, B — V. 
From this we derive the V- and B-band magnitudes, My 
and Mb, respectivel y. Finally, we deriv e the B — I colour 
using the relation of iNatali et al.l (1 19941 ) and with this the 
/-band magnitude, Mi, as well as the V — I colour. Mass 
segregated clusters should show a difference in V — I colour 
in the i nner part with respe ct to the outer part of the cluster. 

(iii) lAllison et al.l (J2009|) suggest a method of detecting 
mass segregation using a minimum spanning tree (MST). 
Their measure, A, shows if a subset of stars is more concen- 
trated compared to a random subset of the same size. It is 
computed with 

A 4 ±J h < 2 » 

where I is the length of the minimum spanning tree of the 
subset, I is the mean length of the MST of random subsets, 
and a is the standard deviation of the distribution of MST 
lengths of random subsets. In a mass segregated cluster the 
most massive stars should show a A well above 1 because 
they are more concentrated than the average subset of ran- 
dom stars. Here we are going to take only stars more massive 
than 1.1 Mq into account, even though the computational 
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expense of this method does barely depend on the total num- 
ber of stars but on the size of the subset. This size, N, was 
practically limited to about 500 in our case, since the com- 
putation time scales with 0(N 2 ) for Prim's algorithm which 
we used (|Primlll957l ). Moreover, we set the n umber of ran- 
dom s ubsets to 50, following the suggestion of lAllison et al.l 

J2009I ) 

(iv) iMaschberger fc Clarke! ([2011]) suggest a different 
method for quantifying mass segregation: by looking at 
the local stellar surface densities of stars. This measure 
has proven to be useful for detecting mass segregation 
in fractal structures which have not merged to a larger 
structure yet, but which ma y already be mass segregated 
IjMaschberger fc Clarke! I2011T ). ft defines mass segregation 
differently than the MST measure: while the MST method 
measures how strongly the massive stars are grouped in re- 
lation to each other, the local surface density method mea- 
sures how strongly the massive stars are grouped in relation 
to all stars. That is, the MST method measures how close 
the massive stars are to each other, whereas the local sur- 
face density method measures how close other stars are to 
the massive stars. We compute th e projected stellar surfac e 
density around each star following ICasertano fc Hut! (jf985), 
i.e. by measuring the radial distance, R, to its 6th nearest 
neighbour in projection and calculating the normalised local 
surface density as 



(3) 



7TjC 2-irnean 

where E me an is the mean local surface density in the cluster. 
In a mass segregated cluster massive stars will have higher 
local densities than the average star. Thus, by normalising £ 
with the mean local surface density we get a dimensionless 
measure which should yield values larger than unity for mass 
segregated stars. For this method we also use only stars more 
massive than 1.1 Mq since the computation of a neighbour 
list for N stars scales with 0(N 2 ) for a brute-force algorithm 
which we use here. To reduce the stochastic scatter, we bin 
the stars in mass bins of 500 stars starting from the most 
massive star. In this way, the A and the local surface density 
measure are better comparable. 



4.2 Substructure 

fn order to detect and quantify substructure and asymmetry 
in our models we test three methods: 

(i) We look at the surface number density profile of mas- 
sive stars. Like for the detection of mass segregation, we 
assume a reasonable completeness level above 1.1 Mq and 
count only stars more massive than that. Inhomogeneities 
will appear as bumps and kinks in this kind of profile. 

(ii) By measuring th e azimuthal density profile (see e.g. 
ICutermuth et al.l 120051) we want to address p ossible asym- 
metries as was done bv lCampbell et alj (|2010r ) for R136. For 
this purpose we count the stars with masses above 1.1 Mq 
within a projected radial distance of 7 pc in 20 azimuthal 
bins of 18 degree each. We quantify the asymmetry by com- 
puting the mean number density of the bins and the stan- 
dard deviation. A comparable measure of asymmetry is then 
given by the normalised standard deviation, i.e. the standard 
deviation divided by the mean. 
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Figure 3. Mass function slopes of the stars with masses above 
1.1 Mq of the models with mass segregation and 100% bina- 
ries among high mass stars. The points show the slopes for 
all stars within a projected radius of 3 pc (thick black) and 
for all stars outside this radius (thin red). The slopes were de- 
termined with the Modified Maximum Likelihood estimator of 
IMaschberger fc Kroupal l|2008r ). The uncertainties were estimated 
using a Monte Carlo approach. For comparison, the green dotted 
line gives the overall mass function slope of 2.3. The red solid line 
is Eq. |4]with coefficients a = —0.09 and b = 2.15. The difference 
in a becomes only evident for strongly mass segregated clusters 
with S > 0.7. 



(iii) ICartwright fc Whitworthl (J200J) suggest a measure 
for determining the degree of substructure, Q. For this 
method, the mean edge length, m, of a minimum spanning 
tree connecting all stars in the cluster has to be measured 
and divided by the mean separation between the cluster 
stars, s. For a homogeneous stellar distribution Cartwright 
& Whitworth find a Q value of about 0.8. A more substruc- 
tured cluster should show a lower value of Q, whereas higher 
values of Q indicate a radial distribution of the stars. 



5 RESULTS 

We test our methods with our calibration models in order to 
see how sensitive the various procedures are in determining 
mass segregation or substructure. In a follow-up investiga- 
tion we aim at applying these methods to A-body models in 
order to see how mass segregation and substructure evolve 
with time in an R136-like configuration. 



5.1 Mass segregation 

5.1.1 Mass function slope 

In Fig. [3] we show the results of the mass function slope de- 
termination for all stars above 1.1 Mq within a projected 
radius of 3 pc and outside this radius. Also shown are the 
uncertainties of these values, which have been estimated us- 
ing 100 Monte-Carlo subsets for each model. While there 
is barely any change in the measured mass function slope 
within 3 pc, at radii larger than 3 pc a changes significantly 
for high degrees of mass segregation, i.e. S > 0.7. At radii 
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Figure 4. The same as Fig.[3]but for the models without binaries. 
The red solid line is Eq. [4] with coefficients a = —0.04 and b = 
2.23. The difference to the models with binaries comes from the 
way the models are set-up: binaries are treated as single particles 
with the sum of the component masses in the set-up process, i.e. 
the difference between the most massive and the least massive 
particle is larger in the case of a high binary fraction and such is 
the degree of mass segregation for any value of S. 



larger than 3 pc the change in a follows a simple relation of 
the form 



a{S) = a(S - l)" 1 + b 



(4) 



with fitted values a = —0.09 and b — 2.15. 

We did the same for the models without binaries 
(Fig. 2| . The change in a with growing S is less pronounced 
compared to the models with 100% binaries among high 
mass stars (a — —0.04 and b — 2.23). This is due to the 
set-up process within McLuster, which first generates the 
binaries, replaces them by centre-of-mass particles with the 
combined mass of the two companion stars, and finally dis- 
tributes those particles within the cluster before they are 
replaced by their constituent stars. With a high binary frac- 
tion the spread in masses between low-mass and high-mass 
particles is higher, and the number of particles is lower dur- 
ing the distribution process. Thus, the final degree of mass 
segregation is higher in the case of high binary fractions. 
For the binary free models, mass segregation becomes only 
significant for models with 5* ^ 0.8. 



5.1.2 Colour gradient 

In Fig. [5] we show the V-I colour profiles for the mass seg- 
regated models. A significant magnitude difference of more 
than 0.1 mag between the inner part of the clusters and the 
outer part is only observable for very high degrees of mass 
segregation (5* > 0.9). 

In Fig. [6] the same is shown for the models without bi- 
naries. Tike for the mass function slope, the effect is less 
pronounced due to the lower effective degree of mass segre- 
gation for models without binaries. In both cases, with and 
without binaries, mass segregation would only be detectable 
for clusters with mass function slopes a ^ 3 outside the core 
using this method on a R136-like young massive cluster. 

It has to be kept in mind, though, that the adopted 
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Figure 5. V -I colour profiles of the m ass segregation models as 
suggested bv lGaburov &: Gielesl I l200a) for the detection of mass 
segregation. Only very high degrees of mass segregation (S ^ 0.9) 
produce gradients larger than 0.1 mag within a radial range of 
about 10 pc. 
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Figure 6. The same as Fig.[5]but for the models without binaries. 
As for the mass function slope (Fig. [4}, the effect of growing S is 
comparable to the case with the high binary fraction. 



upper initial mass lim it in our models is 100 Mq , whereas 
ICrowther et alj (|2010T ) estimate the number of stars exceed- 
ing this limit to be of the order of 10. Such high-mass stars, 
which may even reac h masses of up to 320 M in R136 
IjCrowther et al.ll2010T ), will contribute significantly to the 
blue part of the spectrum, thus will make an observable 
colour gradient more likely if those stars are preferentially 
found near the cluster centre. 



5.1.3 Minimum spanning tree 

In Fig. [7] the min imum spanning tree measure, A, of 
Alliso n et a.1.1 ((2009) is shown for some of the mass segre- 
gated models. The mean MST length was determined for 
all stars in bins of 500 stars and divided by the mean MST 
length of 50 random subsets of 500 cluster stars each. The 
dotted lines give the standard deviation of the random sub- 
sets from this mean value (see eq. [2j. It demonstrates the 
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Figure 7 . Minimum spannin g tree (MST) measure, A, as sug- 
gested by lAllison et al.l II2009T ) . The thicker lines show the nor- 
malised MST lengths for bins of 500 stars versus mean stellar 
mass, <m>. The thinner lines show the standard deviations from 
the mean of 50 random subsets of stars, which is a measure of the 
significance of the detections. All models show a jump at 5Mq 
which corresponds to the mass limit of binaries with similar mass 
companions. Below this threshold the binaries are paired ran- 
domly. The measure shows a slowly but continuously increasing 
degree of mass segregation for all models with S ^ 0.95 and an 
extreme behaviour for the completely mass segregated model. 
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Figure 8. The same as Fig.[7]but for the models without binaries. 
The curves do not show the jump at 5Mq evident in the models 
with the high binary fraction. Only high values of S ^ 0.9 yield 
a significant signal. As in the case with binaries, the difference 
between complete mass segregation and lower values of S is again 
high. 



large significance of the detected mass segregation in all clus- 
ters. 

We see that all curves have a jump at 5 Mq , even the 
5* = 0.0 case, which is due to the binary component pairing 
routine in McLuster. We chose to pair massive O- and B- 
stars with similar mass companions, whereas all stars with 
masses less than 5 Mq are paired randomly. This affects the 
MST length of the massive stars significantly. In Fig.[8]we see 
that this jump disappears for the models without binaries. In 
both cases, most curves show a similar behaviour (A mostly 
between 1 and 3) with a clear trend to higher A for high 
values of S. Only the curve for S = 1.0 significantly stands 
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Figur e 9. Local surfac e density measure, £ (eq. 3; 
iMaschberger &: Clarke 2011). The lines show the normalised, me- 
dian local surface density versus stellar mass. The lines were aver- 
aged with bins of 500 stars. This measure is only weakly affected 
by the binary pairing and shows a smooth transition from unscg- 
regated to completely mass segregated models. 
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Figure 10. The same as Fig. |9]but for the models without bina- 
ries. On average the curves lie at lower values of S and show a 
smaller scatter than for the models with binaries. 

out from the rest. This is due to the fact that in the case 
of S — 1.0 the high-mass stars have the smallest MST that 
can possibly be made. For values of S smaller than 1.0 it 
becomes likely that larger edges are included in the high- 
mass MST such that its length grows rapidly. 

5.1.4 Local surface density 

In Fig. [5]we show the normalised local surface density, E, for 
bins of 500 stars. The curves look similar to the MST curves 
but appear to be less influenced by the ordered binary pair- 
ing. Comparing Fig. [5] with Fig. 1101 which shows the same 
measure but for the clusters without binaries, shows that the 
binaries have indeed only a minor influence on the curves. 
In both cases, the E-values of the highest mass bins grow 
continuously with increasing S. Below ~ 10 Mq the local 
surface density measure suffers from statistical variations, 
though. 

In contrast to the MST measure, the local surface den- 
sity measure does not show a jump between the S = 0.95 
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Figure 11. Summary of the methods for quantifying mass seg- 
regation for the uppermost mass bin (500 most massive stars) 
of all models with mass segregation (solid lines are models with 
binaries, dashed lines are for models without binaries). The A is 
the MST measure and the S is the local surface density measure. 
The ordinate gives the normalised MST length or the normalised 
local surface density, respectively. Note: whereas A = £ = 1 sig- 
nifies no mass segregation of the highest mass stars, A = X is not 
equivalent to £ = X in general. 



model and the completely mass segregated model. This is 
due to the fact that the local surface density of the high- 
mass stars is less influenced by outliers, since it is the av- 
erage of 500 surface density values, whereas the high-mass 
MST length is a sum of 500 edge lengths where one outlier 
can make a large contribution. 



5.1.5 Comparison 

We have shown that the mass function slope method ap- 
pears to be the easiest and most reliable way to detect mass 
segregation for a rich cluster like R136, even though it is not 
very sensitive to low degrees of mass segregation and suffers 
from the arbitrariness of the choice of radius. That is, the 
results of this method depend on the radial range in which 
the mass function is measured and on the underlying mass 
function of the cluster. This complicates the interpretation 
of its results. 

The colour gradient method is also easy to work with 
but is very insensitive in an R136-like configuration. A fur- 
ther test including stars with very high initial mass (Js 
100 Mq) would be very interesting, though, since, due to 
crowding and incompleteness effects, all but the colour gra- 
dient measure are almost impossible to compare with obser- 
vational data of such a cluster. 

The MST measure and the local surface density method 
do work for a R136-like cluster for non-seeing limited data, 
but are computationally much more demanding than the 
other two methods. In contrast to the others, both measures, 
MST and local surface density, allow not only to detect but 
also to quantify mass segregation. That is, their results au- 
tomatically relate the behaviour of the massive stars to the 
other stars. In practice, the quantification is complicated by 
stochastic fluctuations. 

The local surface density measure has the advantage 



w 



100000 

10000 

1000 

100 

10 




0.01 



0.1 



10 



R 



Figure 12. Radial number density profile for stars more massive 
than 1.1 Mq of the clusters with fractal substructure. For lower 
values of D the profiles deviate more strongly from the spherically 
symmetric case (D = 3.0). Only the model with the highest degree 
of substructure (D = 1.6) shows a significant bump at about 0.5 
pc and a kink in the power-law profile. 



that, once the neighbour density of each star is calculated, 
the stochastic fluctuations can be reduced easily by increas- 
ing the bin size which, in addition, increases its sensitiv- 
ity. This is not possible with the MST measure, since new 
minimum spanning trees have to be constructed when the 
sample size is changed. With a sample size of 500 the local 
surface density measure shows a smoother behaviour in the 
very high-mass part, whereas the MST measure is, on aver- 
age, smoother down to lower masses. This may complicate 
the interpretation of the local surface density measure when 
looking at the whole mass spectrum. Moreover, the MST 
measure has the advantage that the significance of its re- 
sults is calculated simultaneously. Therefore, an estimate of 
the significance of the local surface density measure's results, 
like the standard deviation of the MST measure, should be 
constructed in order to make the local surface density mea- 
sure more valuable. 

In Fig. 1111 we compare the MST measure with the local 
surface density method. In this figure we only show the value 
of the uppermost mass bin, and show its dependence on the 
degree of mass segregation, S. Both measures show a steep 
rise for the highest values of S. Moreover, both measures are 
affected by a high binary fraction which may be due to the 
way we construct our binary population, since we pair mas- 
sive stars above 5 Mq preferentially with each other. Never- 
theless, the effective degree of mass segregation should grow 
monotonically with increasing 5*. In this respect, the local 
surface density measure shows a more monotonic behaviour, 
even though it also suffers from statistical fluctuations (pink 
lines). It also gives a value of about 1 for the completely un- 
segregated cluster in the case of a high binary fraction (pink 
solid line). 



5.2 Substructure 

5.2.1 Radial density profile 

In Fig. [12] we show the projected 2d radial number density 
profiles of the clusters with fractal substructure. For the 
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Figure 13. Projected azimuthal number density profile for stars 
more massive than l.lMg of the clusters with fractal substruc- 
ture. The density was measured in 20 bins of 18 deg for all stars 
within a projected radial distance of less than 7 pc from the clus- 
ter centre. The azimuthal variations grow for lower values of D, 
i.e. for higher degrees of substructure. 
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Figure 15. Relative azimuthal variations of the clusters with 
fractal substructure. Shown are three different sets of models each 
with a different random seed. 

The solid line gives a linear fit to all the data (Eq. [5J with a slope 
of -0.46 and an intercept of 1.45. 




Figure 14. Sketch of how the projected azimuthal number den- 
sity profile is determined. The density is measured in 20 bins of 
18 deg for all stars within a projected radius of 7 pc from the 
cluster centre. 



profiles we only took stars more massive than 1.1 Mq and 
put them into 20 logarithmically spaced bins between 0.05 
pc and 20 pc radial distance. A cluster without substructure 
(D = 3.0) is shown for comparison (red solid line). We see 
that for higher degrees of substructure the deviations from 
the spherically symmetric cluster grow. 

For the cluster with D = 1.6, which has the highest de- 
gree of substructure in our sample, the radial profile shows 
a bump and a kink in the slope at about 0.5 pc. Depend- 
ing on the radial range, a power-law fit of the profile would 
yield quite different results for such a cluster or the profile 
could be even interpreted as following a two-part power-law. 
For lower degrees of substructure, the radial density pro- 
file shows much less pronounced deviations making this plot 
rather unfeasible for quantifying the degree of substructure. 



5.2.2 Azimuthal density profile 

In Fig. [13] we show the azimuthal density profiles of the same 
clusters as in Fig. 1121 For this purpose, we used again only 
stars with masses above 1.1 Mq. Moreover, we counted only 
stars within a projected radius of 7 pc and put them into 



20 bins of 18 degree width each (see Fig. ll4|l . The bins were 
then divided by the covered area in pc 2 . Here, the different 
degrees of substructure are more apparent. The spherically 
symmetric model (red solid line) shows only some minor 
statistical fluctuations, whereas the model with the highest 
degree of substructure shows a difference of about a factor 6 
between the bin with the fewest and the bin with the most 
stars. The mean azimuthal density is for all clusters about 
the same (see also Tab. [TJ. 

In Fig. [T^Jthe relative azimuthal variations for all three 
sets of substructured models are shown, that is, the standard 
deviation of azimuthal densities from the mean azimuthal 
density, d£, divided by the mean, E. From the different sets 
of models which are generated with three different random 
seeds we can see that the generation of models with sub- 
structure is a quite stochastic process. But, as expected, the 
relative azimuthal variation grows with decreasing values of 
D. The growth follows a simple relation of the form 



(D) = aD + b, 



(5) 



d£ 

IT 

with a ~ -0.46 and b ~ 1.45. 

The azimuthal density profile appears to be a good mea- 
sure for substructure in star clusters which allows a reason- 
able quantification. 



5.2.3 Q parameter 

In Fig. 1161 we show the Q parameter as defined by 
ICartwright fc Whitworthl (|2004m for our models with fractal 
substructure. The computation of this parameter necessi- 
tates construction of a minimum spanning tree of all cluster 
stars. As stated above, this is computationally very expen- 
sive for a cluster of 170,000 stars when using Prim's algo- 
rithm. Even for a subset of about 15,000 stars with masses 
larger than 1.1 Mq the computation of this parameter takes 
several hours on a regular workstation. One would have to 
switch to a more sophisticated algorithm to use the full set 
of stars which is beyond the scope of this paper. We there- 
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Table 1. Azimuthal number density variations of the three sets of substructured models. £ gives the mean azimuthal number density 
of stars above 1.1 Mq within a projected radial distance of 7 pc, whereas d£ gives the standard deviation from this mean. 



D 



E [pc- 



Sct 1 
dS [pc- 2 ] 



dS/S £ [pcT 



Sct 2 
dS [pc- 2 ] 



dE/S £ [pcr 



Set 3 
dS [pc- 2 ] 



d£/£ 



1.60 


55.2 


33.4 


0.61 


46.5 


34.1 


0.74 


42.4 


32.9 


0.79 


1.70 


66.0 


43.1 


0.65 


32.5 


22.3 


0.68 


49.1 


34.9 


0.71 


1.80 


53.6 


35.9 


0.65 


39.6 


28.4 


0.73 


52.5 


37.7 


0.72 


1.90 


54.4 


29.9 


0.55 


40.6 


20.2 


0.50 


47.2 


25.3 


0.54 


2.00 


41.7 


24.3 


0.58 


38.1 


17.6 


0.47 


42.7 


25.7 


0.61 


2.10 


54.1 


22.2 


0.41 


49.3 


31.6 


0.64 


36.7 


17.4 


0.48 


2.20 


46.1 


18.5 


0.40 


45.3 


22.2 


0.49 


52.8 


20.6 


0.39 


2.30 


46.7 


18.3 


0.39 


51.0 


17.5 


0.36 


57.4 


22.3 


0.39 


2.40 


45.5 


18.6 


0.41 


47.8 


19.2 


0.40 


47.7 


12.3 


0.25 


2.50 


52.4 


14.3 


0.27 


44.4 


14.4 


0.33 


55.4 


15.6 


0.29 


2.60 


49.7 


11.6 


0.23 


48.4 


10.1 


0.21 


40.7 


11.4 


0.28 


2.70 


43.9 


10.3 


0.23 


46.9 


8.70 


0.19 


43.8 


11.3 


0.26 


2.80 


44.2 


8.02 


0.18 


46.4 


14.7 


0.32 


44.4 


8.19 


0.19 


2.90 


46.0 


6.16 


0.13 


47.5 


7.25 


0.15 


49.8 


6.51 


0.13 


3.00 


47.7 


3.20 


0.07 


47.0 


3.30 


0.07 


47.1 


3.98 


0.08 



o 



1.4 
1.2 

1 
0.8 
0.6 
0.4 
0.2 h 





i 1 1 1 1 1 1 r 



%fH0T- 




+ + 



m>1.1 M (Set1) + 

m>1.1 M (Set2) X 

m>1.1 M„(Set3) )K 

m>5.0M o (Set1) + 



1.6 1.6 



2.2 2.4 
D 



2.6 2.8 



Figure 16. Cartwright &c Whitworth's Q parameter of the clus- 
ters with fractal substructure. Computation of this parameter is 
computationally too demanding when all stars are taken into ac- 
count and Prim's algorithm is used for the computation of the 
MST length (see text). But when calculating it for only a subset 
of the most massive stars, the absolute value of the Q parameter 
seems to depend on the number of stars in the subset, which is 
most probably due to the binary population. The lower solid line 
gives a linear fit to the Q parameters for the subsets containing 
all stars above 5 Mq (Eq. [7)l with a slope of 0.16 and an intercept 
of 0.44. In contrast, the upper solid line is a fit to the subsets 
containing all stars with masses larger than 1.1 Mq (about 15000 
stars) from all three sets of models. The slope here is 0.18 and 
the intercept 0.75. 



fore compute the Q parameter for different sizes of subsets 
with different low-mass cut-offs. 

We compute Q by measuring the mean separation be- 
tween the cluster stars, s, and the mean edge length, rn, of 
a minimum spanning tree connecting all stars in the subset, 
such that 

Q=^- (6) 

s 

Cartwright & Whitworth find values between 0.45 for highly 
fractal configurations (fractal dimension D — 1.5) and 0.8 



for homogeneous configurations (D — 3.0). Values above 
0.8 and up to about 1.5 they find for models with a radial 
density gradient {pzd oc r~ s ). But, in contrast to their in- 
vestigation, the models in our sample have a radial density 
gradient {p2d oc R~ 2 ) and at the same time fractal sub- 
structure (1.6 ^ D ^ 3.0). Therefore, we expect the value of 
Q of our models to be of order of the values Cartwright & 
Whitworth find for the radial clusters, but also to show some 
variability around this value depending on the degree of sub- 
structure. Moreover, our models contain binaries which are 
paired non-randomly, i.e. the massive stars above 5 Mq are 
paired with similar mass companions (see Sec. [3} which fur- 
ther complicates the interpretation of Q. 

Interestingly, as can be seen in Fig. 1161 the Q values 
of our models depend on the size of the subset. For a small 
subset containing only the most massive stars with masses 
larger than 5 Mq (about 2500 stars) the values almost agree 
with the findings of Cartwright & Whitworth for purely frac- 
tal clusters but are shifted about 0.1 upwards. The more 
stars are taken into account the more the absolute values 
of the Q parameter are shifted to larger values. This is due 
to the binaries in our clusters: when taking all stars above 
5Mq we only have massive binaries in our sample, as the 
massive stars are paired with each other. When going to a 
lower mass limit we add stars to the sample that need not 
be in binaries, i.e. whose binary companions need not to be 
within the sample. This increases the mean edge length rn, 
whereas it barely affects the mean separation s, such that Q 
becomes larger. 

We fit a linear function to the values of the form 



Q(D) = aD + b, 



(7) 



where a = 0.16 and b — 0.44 for the small subsets with 
m > 5 Mo, whereas a — 0.17 and b = 0.75 for the large 
subsets with m > 1.1 Mq (about 15000 stars) from all three 
sets of models. Subsets with sizes in between these two yield 
intermediate values (e.g. a — 0.21 and b — 0.62 for a subset 
with m > 1.9 Mq, i.e. about 7500 stars). 

To which exact values Q will converge when more stars 
are taken into account cannot be checked easily but the 
growth in Q with increasing sample size falls off strongly 
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such that the overall values should be close to the values 
of our largest samples. Anyway, the question remains how 
useful such a quantification of substructure is. To follow the 
evolution of substructure with time in a single JV-body com- 
putation seems possible with this method. But it is time 
consuming, unless a highly sophisticated algorithm is used. 



6 SUMMARY AND CONCLUSIONS 

Here we introduce for the first time the code McLuster, 
which is a publicly-available tool for initialising star cluster 
models and binary-rich stellar populations. 

Moreover, we have tested and calibrated several meth- 
ods for detecting and quantifying mass segregation and sub- 
structure in non-seeing limited star cluster data. We applied 
these methods to models of the young massive cluster R136 
which is the only starburst site in the Local Group. Access- 
ing this cluster and measuring its degree of mass segregation 
and substructure at its current age of about 3 Myr is of im- 
portance because it is the only object which can be resolved 
into single stars and at the same time can give insight to star 
formation conditions in starburst sites. Moreover, due to its 
size and mass we can assume that the present-day condi- 
tions in R136 can still yield insights to the conditions 3 Myr 
ago, as the dynamical time of this configuration is compa- 
rably l ong (about 0.3 Myr), i.e. the cluster is dy namically 
young jPortegies Zwart. McMillan fc Gielesl l201ol l. Finally, 
its moderate mass of about 10 5 M© makes it accessible by 
means of iV-body investigations. But for understanding the 
development of these two quantities, mass segregation and 
substructure, easily quantifiable and little time consuming 
methods have to be found. 

We compare 4 different methods for quantifying mass 
segregation and 3 methods for quantifying substructure from 
the literature (Sec. 2J. The former we quantify by com- 
paring the m ass function slope of ma ssive stars, the radial 
colour profile lIGaburov fc Giclcs 2008), the minimum span- 
ning tree measure (I Allison et al.l 2009. 20 10]) and the lo - 
cal surface density measure ( Maschberger fc Clarkej|201ll ). 
We quantify substructure by looking at the projected ra- 
dial nu mber density profile, the proj ected azimuthal density 
profile JGutermuth et al|2 005. 20 08|) and by calculating the 
Q parameter ( Cartwright fc Whitwortbll2004l ) . For this pur- 
pose we set up star cluster models with different degrees of 
mass segregation and substructure using the new publicly 
available code McLuster (Sec. [3j . 

We find that the methods for detecting mass segregation 
often yield ambiguous results (Sec.[5)|. From the four meth- 
ods we compare, the mass function slope seems to be the 
simplest and most reliable measure for detecting mass seg- 
regation, even though it suffers from a somewhat arbitrary 
choice of radius, and significant detections of mass segrega- 
tion are only possible for high degrees of mass segregation. 
The radial colour profile only yields significant results for the 
highest (and rather unrealistic) degrees of mass segregation 
for a configuration like R136. Both measures cannot handle 
substructure, i.e. clusters which are not radially symmetric. 
In contrast to that, the minimum spanning tree measure and 
the local surface density measure can handle substructured 
clusters (JMaschberger fc Clarkdl201ll ). but their results are 
often ambiguous, and computationally they are much more 



demanding. Moreover, the minimum spanning tree method 
is strongly influenced by a high binary fraction. On the other 
hand, it has the unique advantage that the significance of 
its results is readily given. The local surface density measure 
is less influenced by binaries, and with careful calibration it 
can be a very sensitive method for detecting mass segrega- 
tion. We recommend, though, that a measure of its signifi- 
cance similar to the standard deviation of the MST measure 
should be constructed. 

For quantifying substructure we are left with the pro- 
jected azimuthal density profile since the projected radial 
density profile only shows significant deviations from the 
spherical symmetric case for extremely substructured con- 
figurations. Such a cluster with a high degree of substruc- 
ture can show strong bumps and kinks in its projected ra- 
dial profile but those are difficult to quantify. The projected 
azimuthal density profile is a reliable measure of substruc- 
ture (Fig. I15p . We suggest to compute the mean projected 
azimuthal density and the normalised standard deviation 
from this mean to get a useful measure. The Q parameter 
is also sensitive to substructure but is computationally very 
demanding, when using a standard algorithm like Prim's al- 
gorithm, as a minimum spanning tree for all cluster stars has 
to be constructed. If only a subset of stars is taken account 
for its computation, its absolute value shows a dependence 
on the number of stars in this subset which is due to the bi- 
nary population we adopted. Thus, the Q parameter is only 
of limited use in such configurations, unless a much more 
sop histicated algorithm like Kruskal's algori thm is used (see 
e.g. lLomax, Whitworth fc Cartwrightlboill ). 

Finally, we have to conclude that most of the methods 
presented in this work will likely yield ambiguous results 
when applied to observations of young massive clusters, due 
to crowd ing and incompleteness effects, as has similarly been 
found bv lAscenso. Alves fc Lagol (J2009J) . Also the distance of 
such clusters and heavy, variable extinction will add further 
difficulties. But even for non-seeing limited data of (young) 
massive clusters, like JV-body data, some of the methods will 
face great computational difficulties due to the large number 
of stars involved. 
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APPENDIX A: MCLUSTER MANUAL 

The tool McLuster is an open source programme that 
can be used to either set up initial conditions for N- 
body computations or, alternatively, to generate artificial 
star clusters for direct investigation. There are two dif- 
ferent versions of the code, one basic version for generat- 
ing all kinds of unevolved clusters (in the following called 
mcluster) and one for setting up evolved stellar popula- 
tions at a given age. The former is completely contained in 
the C file main. c. The latter (here dubbed as mcluster_sse) 
is more complex and requires additional FORTRAN rou- 
tin es, namely the Sing l e-Star Evolution (SSE) routines 
by iHurlev. Pols fc Tout! (J2000h which are provided with 
the McLuster code. For a quick introduction read the 
README file which is also provided with the code. For tech- 
nical details on how to generate initial con ditions f o r star 
cluster in general we would like to refer to iKroupal 12008) 
and referenced literature therein. 



Al Compilation 

After extracting the archive which can be obtained from 
the given web addresqj, the basic version mcluster can 
be compiled on a UNIX system from the command line with 

> cc -lm -o mcluster main.c 



www . astro . uni-bonn . de/~akuepper/mcluster/mcluster . html 
or www . astro . uni-bonn . de/~webaiub/english/downloads . php 



where cc may be replaced by the C-compiler avail- 
able on your computer. It can also be compiled by using 
the Makefile, i.e. 

> make mcluster 

For the more complex version mcluster_sse the Makefile 
has to be used. Type 

> make mcluster_sse 

after which the programme should compile, generating 
an executable named mcluster_sse. 

• Note that, when using the Makefile, you may have to 
change the C- and/or FORTRAN-compiler entry as well as 
the path of your compiler. 

• In case you want to apply any change to the code make 
sure that you first delete all object files by typing 

> make clean 

before re-compiling the code. 



A2 Input 

There are two ways of choosing the desired cluster parame- 
ters. One is to set the parameters manually within main.c, 
within the upper part of the code right at the beginning of 
the main routine where all variables are declared. Note that, 
after changing the value of a variable, you have to compile 
the code again. The more convenient way is therefore to 
pass arguments to the code at the time of execution via 
the command line (see Tab. [l] for an overview of available 
options). Type 

> mcluster -h 



> mcluster_sse -h 

respectively, to get a quick help on the available pa- 
rameters and their usage. 

• In case you have not specified a certain parameter, the 
default value as set within the main routine is used. 

• Not all parameters can be set via the command line, 
some have to be changed within the code. 

• All command line arguments are the same for mcluster 
and mcluster _sse, except for the age parameter (-e) which 
is mentioned in more detail below. 



A3 Density profile 

McLuster can generate star clusters with various radial 
density profiles (option -P). In all cases the total mass of 
the cluster has to be chosen separately (options -M or -N). 
The available profiles are: 

(i) The simplest option is the analytical IPlummerl (|1911h 
profile (option -P0) which is the simplest (two-parameter 
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Table 1. Overview of available command line options in McLuster with brief descriptions. For details on the available choices see the 
corresponding paragraphs. Also given are the possible ranges and the default values (which can be permanently changed within main.c). 



Option 



Range 



Default 



Meaning 



Number of cluster stars (specify either N or M) 

Mass of the cluster (Mq; specify either N or M) 

Density profile (Plummcr/King/Subr/EFF/homogeneous sphere) 

Wo parameter for the King model (only P = 1, ignored if P ^ 1) 

Half-mass radius in parsec (ignored for P = 3) 

Scale radius of the EFF template in parsec (only P = 3, ignored if P ^ 3) 

Cut-off radius of the EFF template in parsec (only P = 3, ignored if P ^ 3) 

Power-law slope of the EFF template (only P = 3, ignored if P ^ 3) 

Degree of mass segregation (0.0 means no segregation; S < 1.0 for P = 2) 

Fractal dimension (3.0 means no fractality) 

NBODY4/6 computation time in N-body units 

Virial ratio (Q = 0.5 for virial equilibrium) 

Output type (Nbody6/Nbody4/ASCII table) 

NBODY4/6 adjustment time in N-body units (e.g. Hcggic & Hut 2003) 

NBODY4/6 output time in N-body units 

Use GPU with Nbody6 (no/yes) 

Output name of the cluster model 

IMF (no mass function/Kroupa IMF/user defined) 

IMF slope for a user defined IMF, may be used multiple times, from low to high mass 

IMF mass limit ( Mq) for a user defined IMF, may be used multiple times, from 

low to high mass 

Number of binary systems (specify either B or b) 

Binary fraction (specify either B or 6) 

Binary pairing (random/ordered for M > 5Mq) 

Seed for randomisation, set for randomisation by local time 

Tidal field (no tidal field/near-field approximation/point-mass potential/ 

Milky- Way potential) 

Epoch for stellar evolution in Myr (only available in mcluster_sse) 

Metallicity (Z = 0.02 for solar metallicity) 

Galactocentric radius vector in parsec (use 3 times for x-, y- and z-coordinate) 

Cluster velocity vector in km/s (use 3 times for x-, y- and z-coordinate) 

Output units (TV-body units/astrophysical units) 

Display help 



-N 


< N < Nraax 





-M 


M > 


1000 


-P 


0/1/2/3/-1 





-W 


1.0-12.0 


— 


-F, 


R> 


0.8 


-r 


< r < c 


— 


-c 


c > r 


— 


-g 


g> 1.5 


— 


-S 


0.0-1.0 


0.0 


-D 


1.6-3.0 


3.0 


-T 


T> 


100.0 


-Q 


Q^0 


0.5 


-c 


0/1/3 


3 


-A 


A > 


2.0 


-O 


O > 


2.0 


-G 


0/1 





-o 


— 


test 


-f 


0/1/2 


1 


-a 


a < 


— 


-m 


m > 


— 


-B 


O-N/2 





-b 


0.0-1.0 


0.0 


-P 


0/1 


1 


-s 


s ^ 





-t 


0/1/2/3 


3 


-e 


e ^ 





-Z 


0.0001-0.03 


0.02 


-X 


X ^ 


8500/0/0 


-V 


V >0 


0/220/0 


-u 


0/1 


1 


hi-! 


— 


— 



only) stationary solution to the collisionless Boltzmann 
equation. For this profile only the half-mass radius has to be 
specified additionally (option -R, in parsec). The Plummer 
profile is in principle infinitely extended but gets automati- 
cally truncated at the theoretical tidal radius of the cluster 
in case a tidal field has been specified (see below). 

(ii) A more sophistic ated s e t of m odels is given by the 
distribution function of iKing (|196q ). For this profile (op- 
tion -Pi) the half- mass radius and the dimensionless value 
of Wo which specifies the model concentration (option -W 
ranging from 2.0, i.e. low concentration, to 12.0, i.e. high 
concentration) has to be specified. The underlying routine 
for generating the density distribution from the distribu- 
tion function is based o n kingO by Douglas C. Heggie (e.g. 
iHeggie fc Aarsethlll992i ). Note that the final density distri- 
bution gets scaled to exactly match the desired half-mass ra- 
dius; the radius at which the density becomes zero does not 
necessarily match the theoretical tidal radius in this case. 



(iii) ISubr, Kroupa. fc Baumgardtl (|2008l ) give a density 
profile which can be chosen to be mass segregated. For this 
density profile (option -P2) the half-mass radius and an 
additional mass-segregation parameter (option -S, ranging 
from 0.0-0.99, but S < 0.5 for reasonable models in virial 
equilibrium) have to be chosen. For -S0.0 the Subr profile 
is equal to a Plummer profile. McLuster uses a slightly 



modified version of the Plumix routine by Ladislav Subr to 
set up this kind of profile. 

(iv) Young clusters in the Large Magellanic Cloud were 
found to follow a two-dimensional density profile consisting 
of a co re and a power-law ta i l with out visible tidal trun- 
cation. lElson. Fall fc Freeman! 1)19871 ) give a simple analyti- 
cal formula for such profiles which can be deprojected with 
McLuster and used to set up 3d star cluster models (option 
-P3). Since those models are in principle infinitely extended, 
the so-called EFF models do not get scaled to a certain half- 
mass radius, but rather require specification of a cut-off ra- 
dius to which the profile should extend (option -c, in par- 
sec) . The central density then gets calculated automatically 
using the specified cluster mass (see below) . In addition, the 
radius of the two-dimensional core (option -r, in parsec) 
and the slope of the power-law part of the profile (option 
-g) have to be chosen. Values for g larger than 1.5 usually 
yield stable solutions, observational values lie at about 2.0 
for y oung star clusters in the LMC IjElson. Fall fe Freeman! 
Il987l ). The velocity field of this fa mily of profiles is generated 
using the algorithm described in Kroupa! l|2008h . 

(v) A further possibility to set up the density distribu- 
tion is given by option -P-l. In this case the final cluster 
has no density gradient, but consists of stars which are ho- 
mogeneously distributed within a sphere. The size of the 
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sphere is specified by choosing the half-mass radius (thus, 
the limiting radius of the sphere will be a factor 2 1 ' 3 larger). 
This option is especially useful in case of fractal initial con- 
ditions (see below). The velocities of the stars are isotropic 
and drawn from a Gaussian distribution. 

• The exact matching of the actual half-mass radius as 
resulting from the discretized model to the specified value 
may be switched off within the main routine, but is not rec- 
ommended. Set match = and compile the code again. 



A4 Tidal field 

As mentioned above, the choice of the tidal field and of the 
cluster orbit may influence the extent of the density profile. 
McLuster offers different kinds of tidal fields which can be 
specified with the option -t. In addition it may be necessary 
to specify the galactocentric radius vector and the orbital 
velocity vector. This can be done by using the option -X 
(in parsec) and -V (in km/s), respectively, three times for 
the x-, y- and z-component (within a Cartesian coordinate 
system where the galactic disk would lie in the x-y-plane, if 
applicable). As an example, -X8500.0 -X0.0 -X0.0 -V0.0 
-V220 . -V0 . would give the motion of the Local Standard 
of Rest. This option is especially useful when generating 
input for JV-body computations, since the input parameter 
file is automatically adjusted accordingly. 

(i) For a cluster in isolation choose -to. No truncation is 
applied to profiles like the Plummer profile in this case. 

(ii) A linearized tida l field, as described in 
iFukushige fc Heggiej (|2000r ) can be chosen with -tl. If 
you have selected to generate input files for Nbody6 (see 
below) then the values for the Local Standard of Rest are 
used. In all other cases the galactocentric distance has to be 
specified additionally. Therefore use the option -X and set 
all but one component to zero, e.g. -X6000.0 -X0.0 -X0.0 
for an orbit at 6 kpc. No orbital velocity vector has to be 
specified as the linearized tidal field mimics a circular orbit. 

(iii) If you choose -t2 then you get a point-mass galaxy 
for which you can specify any galactocentric radius and or- 
bital velocity (options -X and -V). The mass of the galaxy 
is set within the header of the main routine. By default, 
Mlpointmass is set such that you get an orbital velocity of 
220 km/s at a galactocentric radius of 8.5 kpc. 

(iv) For a more realistic, Milky- Way potential you can 
choose option -t3. This potential consists of a ce ntral point 
mass /bulge, modelled as a Hernquist potential (|Hernquistl 
|1990T ). a Miyamoto disk (|Mivamoto fc Nagail Il975h and a 
logarith mic (phantom dark-ma tter) halo, with values as 
given in I Allen fc Santillanl |l99ll 'l. If you set up initial condi- 
tion for Nbody6 then the logarithmic halo will be adjusted 
such that the circular velocity, VCIRC, at some radius, RCIRC, 
has a specific value. The default is 220 km/s and 8.5 kpc, re- 
spectively, which may be changed in the header of the main 
routine. The other parameters of this potential may also be 
changed there. 



A5 Cluster mass and stellar mass function 

You can either fix the total number of stars in your cluster 
(option -N) or you can set a desired total mass (option -M, 



in solar units). In the latter case, McLuster draws stars 
from the selected mass function until the desired mass is 
exceeded. The mass function of stars in the cluster can be 
defined to be one of the following three kinds (option -f ). 

(i) All stars can have the same mass (option -fO). The 
mass of each star is by default assumed to be 1 M© , 
which may be changed within the main routine (parameter 
single_mass). 

(ii) The canonical iKroupal (|200lh initial mass function 
can be used with -f 1. This IMF has a slope of ai = —1.3 
for stellar masses m = 0.08 — 0.5 Mq, and the Salpeter slope 
a 2 = -2.3 for m > 0.5 M . The lower and upper IMF limit, 
mlow and mup, are by default 0.08 Mq and 100 Mq, respec- 
tively, but these values may be changed in the main routine. 

(iii) More sophisticated, multi-power-law IMFs can be set 
up with option -f 2. Therefore, McLuster uses the routine 
mufu by Ladislav Subr. This routine allows to define several 
mass limits and corresponding mass-function slopes between 
these limits. The limits and the slopes can be passed to 
McLuster with the options -m and -a, respectively. These 
options have to be used multiple times, where one more limit 
has to be passed to McLuster than number of slopes. For 
example, the Kroupa IMF with a steeper slope of 03 = —2.7 
for stars more massive than 5 Mq up to a maximum mass of 
80 M Q would be -f2 -m0.08 -a-1.3 -m0.5 -a-2.3 -m5.0 
-a-2.7 -m80.0. 

• The maximum stellar mass - cluster mass relation 
found by IWeidner fc Kroupal i|2006T ) can be used to auto- 
matically cut off the IMF at the corresponding upper mass 
limit. This routine is switched off by default but may be acti- 
vated by setting the weidner parameter in the main routine 
to 1. 

• In McLuster there is a maximally allowed stellar mass 
limit defined through the parameter upper_IMF_limit. This 
parameter is set to 100 Mq since Nbody6, i.e. the stellar 
evolution routine SSE within Nbody6, does not allow higher 
masses. In case you need higher stellar masses anyway, set 
this parameter within the main routine to the desired value 
but keep in mind that stars with mass larger than 100 Mq 
are evolved as 100 Mq stars. 

A6 Mass segregation 

With McLuster it is possible to apply any degree 
of primordial mass segregation to all available density 
profiles, not only to the Subr profile as already men- 
tioned above. Therefore, the method as described in 
iBaumgardt, De Marchi fc Kroupal l|2008h is used. The ad- 
vantage of this method is that the chosen density profile is 
not changed with increasing degree of mass segregation as 
is the case for the Subr profile. 

In short, it works as follows: for a cluster of N stars 
McLuster first draws the stellar masses from the selected 
IMF (see above) and then creates N' = N{m)/mi ow orbits, 
where (to) is the mean stellar mass and mi ow the lowest 
stellar mass in the cluster. These orbits get ordered by their 
specific energy, from low energy to high energy orbits. Then 
the stellar masses are also ordered and the cumulative mass 
function, M cum (i) — X/ = i M(j) is evaluated from this mass 
array, whereby i = 1 . . . N. After dividing M cum (i) by the 
total cluster mass, the function is normalised such that it 
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runs from to 1. Finally, for any star an orbit from the list 
of energy-ordered orbits is chosen from the orbits between 
N'M CU m{i - 1) and N'M cum (i). 

If the masses in the mass array are perfectly ordered 
from highest to lowest then this will yield a completely mass 
segregated cluster. That is, the highest mass star is on the 
lowest energy orbit, and the lowest mass star is on the high- 
est energy orbit. 

Intermediate degrees of mass segregation can be 
achieved by non-perfect ordering. In McLuster this is re- 
alised as follows: first, all N stellar masses are ordered from 
highest to lowest. Then, beginning with the highest mass, 
the masses are written to a new array, where the i-th mas- 
sive star is written to the j-ih empty slot counting from 
to N — i. j is generated using 



j = (N-i)(l-X 1 - 3 ) 



(Al) 



where X is a random number between and 1, and S is 
the mass segregation parameter. When S is zero, j can have 
all values from to N — i and we end up with a random 
distribution. But when S is 1, then j is always zero and we 
reproduce the perfectly ordered array we started with. That 
is, because every star i, beginning with the most massive 
one, gets written to the next empty slot which is slot i. By 
choosing S values between and 1 we can get intermedi- 
ate degrees of partial mass segregation (option -S). Unlike 
for the Subr profile, S can explicitly be chosen to be 1.0. 
Moreover, for all values of S we get clusters in virial equilib- 
rium (if not explicitly specified differently) with the desired 
(mass) density profile. 



A7 Fractality 

Star clusters are not born in perfect symmetry. They are 
rather formed in collap s ing, fractal molecula r clouds (see e.g. 
iGutermuth et alj|2008i : iKonvves et alj|20ld ). With McLus- 
ter you can set up two kinds of fractal initial conditions. 
First, you can set up a fractal distribution of stars within a 
sphere of constant average den sity, similarly as described in 
iGoodwin fc Whitworthl (2004). Secondly, you can add frac- 
tal substructure to any of the above given density profiles. 

(i) When you choose a homogeneous density profile (op- 
tion -P-l) the cluster stars get distributed within a sphere 
as follows. The first star (a so-called parent) is placed in the 
middle of a box of size 2 (arbitrary units), then this box 
is split into 8 sub-boxes of half the initial box size. In the 
centre of each sub-box a further star is placed (a so-called 
child), whereupon each sub-box is split up into 8 smaller 
pieces, such that each child now becomes a parent on its 
own. By applying a small random offset to each star from 
its sub-box centre, we make sure that the final cluster does 
not look too grid-like. This is repeated until we have gener- 
ated 128.0 x 8 log(Ar)/ log(s) stars or until the total number of 
stars would exceed this number with the next generation of 
children. From these stars we randomly draw N stars with 
radial distances of less than unity from the centre of the 
initial box. We end up with a homogeneous sphere of stars. 

Now, if not every sub-box gets a new star, and only those 
sub-boxes get sub-divided which have a star in their cen- 
tre, then the final distribution of stars becomes fractal. The 
probability that a sub-box gets a star can be expressed as 



2^ D ~ 3 \ where D is the fractal dimension (option -D). If D is 
chosen to be 3.0 then we get no fractality since the probabil- 
ity is unity, i.e. every parent gets 8 children. If it is, e.g., 2.0 
then only every second sub-box gets a star, or, on average, 
every parent has 4 children. 

Corresponding stellar velocities are drawn from a Gaus- 
sian distribution, and re-scaled such that all children of one 
parent are in virial equilibrium and the total mean velocity 
in one sub-group is unity (arbitrary units). In addition, each 
child gets the velocity of its parent. In a later step, McLus- 
ter re-scales the phase-space coordinates of the stars such 
that the cluster is in virial equilibrium (if not specified dif- 
ferently). In this way, we get a fractal structure consisting 
of coherently moving, gravitationally bound substructures. 

(ii) Alternatively, we can choose to set up fractal clusters 
which follow a given density profile like, e.g., the Plummer 
profile or the King profile, but which show fractal substruc- 
ture. This is realised by first generating a sphere of stars of 
radius unity with the above procedure. But now this distri- 
bution of stars gets folded with the chosen density profile. 
Therefore the radius of each star first gets re-scaled by its 
absolute value to the power of three. In this way, we get N 
stars with radii distributed homogeneously between and 
1, but which show substructure in 6D phase space. These 
radii are used as seeds to compute a corresponding radius 
within the specified density profile. In a last step, the space 
coordinates of each star get scaled by this newly generated 
radius. In this way the fractal distribution is conserved but 
folded with the specific density profile. Moreover, the veloc- 
ity of each star is scaled to the expected velocity of a star 
at the given radius within the specified density profile. 

Note that the method described here is an ad hoc in- 
troduction of substructure which, like all other methods for 
generating substructured models, does not rely on any phys- 
ical motivation. In this way, the generated clusters can have 
any degree of substructure and a smooth transition from 
spherical symmetry to substructuredness can be realised. 
This may be useful in some applications but is not meant 
to accurately reproduce observational substructure. In fact, 
due to the radial re-scaling the substructure gets stretched 
which produces long filaments in the final clusters. These 
filaments may be compared to the filamentary structure of 
molecular gas in star forming regions. 

• In order to guarantee a minimum of spherical symmetry, 
you can tell McLuster to give 8 children to the first "ur- 
parent" . In this way, you avoid too large asymmetries. This 
will lead to less differing results between initial conditions 
generated with different random seeds, but does on the other 
hand not yield perfectly fractal clusters. This tweak can be 
switched off within the header of the main routine. Set the 
symmetry parameter to and re-compile. 

• Once in a while McLuster gets stuck in the fractality 
sub-routine. In this case a restart with a different random 
seed should help. 



A8 Binaries 

After the stellar masses got drawn from an initial mass func- 
tion, you can choose to let McLuster set up binary systems. 
You can specify either the desired number of binary systems 
(option -B, values from to N/2) or you can specify a frac- 
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tion of stars which should be in binary systems (option -b, 
values from to 1). Thus, from all N stars, N x 6/2 or B 
binaries are formed, respectively. The binaries are then re- 
placed by a centre-of-mass (CoM) particle for the rest of the 
procedure. Only in the very end, after the density profile 
has been established and the velocities of the cluster mem- 
bers have been scaled appropriately, the CoM particles get 
replaced by their two constituent stars. The binary orbital 
plane is oriented randomly and the orbital phase is also cho- 
sen randomly. The internal properties of the binaries can 
be generated according to the following semi-major axis dis- 
tributions which can be selected in the header of the main 
routine (parameter adis). 

(i) A flat distribution of semi-major axes can be speci- 
fied with adis = 0. In addition, you have to choose a min- 
imum and a maximum semi-major axis (parameters amin 
and amax). 

(ii) If adis is set to 1 (default) then the semi-major 
axis of each bi nary is comput ed from a period which was 
drawn from th e lKroupal (|1995aT ) period distribution (see also 
I Kroupal 120081 ). This is the preferred distribution function 
as it unifies the observed Galactic-field and the pre-main- 
sequence population. 

( iii) If you want the semi-maj or axes to be generated from 
the lDuquennov fc Mayor] (| 199 11 ) Galactic-field period distri- 
bution then you have to set adis = 2. 

• The pairing of primary and secondary components of 
the binaries can be chosen to be either random or ordered. 
In the latter case (option -pi) the stellar masses above a 
certain threshold (parameter msort in the main routine) get 
ordered from most to least massive. The rest of the stars are 
put in random order onto the list below the last star with 
mass above msort. The pairing of binaries now starts with 
the most massive star which gets paired together with the 
second-most massive, then the third with the forth, and so 
on down the list. Below msort pairing is random, which is 
consi stent with the binary-star observational data (|Kroupal 
[2003). The default value of ms ort is 5Mq, in ro ugh agree- 
ment with recent findings fe.g. iKobulnickv fc Frverl l2007V 

• Eccentricities, e, are drawn from a thermal eccentricity 
distribu tion, i.e. f(e) = 2e (e.g. iDuquennov fc Mavorlll99ll ; 
see also lKroupall200S| y 

• To correct for the fact that short-perio d binaries in th e 
Milky Way do not show high eccentricities (Mathicu 1994), 
iKroupal |l995bj) introduces an analytical correction for such 
systems, which is attributed to pre-main sequence eigenevo- 
lution between the constituent stars. This correction can be 
switched off by setting the parameter eigen = in the main 
routine. 

A discus sion of binary systems and their formation can be 
found in lKroupal (|2009l ). 



A9 Stellar evolution 

McLuster conta ins the SSE routines by 
Irlurlev. Pols fc Toutl (120001 ) which are also used in, e.g., 
Nbody6 for stellar evolution. If you just want to generate 
star clusters consisting of zero-age main sequence (ZAMS) 
stars then you only need the basic version mcluster. But 
if you want to set up a cluster with an evolved stellar 



population then you have to use mcluster_sse which makes 
use of those SSE routines. In this case you have to specify 
an age for the cluster population (option -e, in Myr). The 
evolution of the stars is done in the very beginning of the 
programme. When the stars are drawn from the IMF they 
get immediately evolved to the desired age. The masses of 
the evolved stars are then summed up and additional stars 
are generated until the desired cluster mass is exceeded 
or the desired number of stars is reached. The stellar 
parameters derived from SSE for each star are stored in an 
additional file, which also has to be passed to Nbody6 (see 
below). 

• The internal parameters of SSE can be changed within 
the header of the main routine (not recommended). 

• In case a star becomes a neutron star or a black hole 
SSE assigns a kick velocity to the remnant (if not specified 
differently). The kick velocity can be used to remove the 
remnant from the cluster. Therefore the present-day escape 
velocity of the cluster at its half-mass radius is calculated 
and if the kick velocity exceeds this velocity it gets removed. 
If you want to keep all compact remnants set the parameter 
remnant = in the main routine. 

• The metallicity, Z, can be set with option -Z. Alterna- 
tively, you can specify the metallicity as [Fe/H] within the 
main routine, the cor responding Z v alue i s computed using 
the relation given in iBertelli et al.l (119941 ) . Make sure that 
in this case you set Z — beforehand. 

• In case you are generating binaries and have selected 
ordered pairing for stars above a certain mass, msort (see 
above), then the ZAMS mass is used to decide whether a 
star is paired randomly to another star or not. 

• The components of binaries are independently evolved 
as single stars with SSE. For a more advanced treatment of 
st ars in binaries, the Binary-S tar Evolution (BSE) routines 
bv lHurlev. Pols fc Toutl (|2002| ) are also included in McLus- 
ter. Therefore, at the very end of the cluster generation 
procedure, when the binary CoM particles are replaced by 
their constituents and the orbital elements are generated, 
the masses of the components are reset to their ZAMS mass. 
Then the two stellar masses, a semi-major axis and an eccen- 
tricity are passed to BSE which finally returns correspond- 
ing evolved values. This feature is switched off by default 
but may be activated by setting the parameter BSE to 1 in 
the main routine. 



A10 Output 

Up to no w McLuster can generate input for Nbody6 (op- 
tion -CO. lAarsethll2003l ) and Nbody4 (option -CI), or it can 
write an ASCII table of stars and their properties (option 
-C3). 

(i) In the first and second case, there will be two out- 
put files which can be named with option -o. For exam- 
ple, -o mycluster will yield the files mycluster .input, 
containing all the input parameters for the run, and 
mycluster .fort . 10, containing the masses, positions and 
velocities. Note that the latter has to be renamed to fort . 10 
at the time of execution in order to be recognised by 
Nbody4/6. When using mcluster_sse there will be another 
file named mycluster .fort . 12. This file also has to be re- 
named within the directory of the run to fort .12. The name 
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mycluster is just added to the file names for convenience. 
Thus, a directory for an Nbody4/6 run should contain: 

(a) mycluster . input, 

(b) fort. 10, 

(c) fort. 12. 

The run is then started with the usual command, i.e., 

> / . . . /nbody6 < mycluster . input 

where / . . . / should be replaced by the path to your 
Nbody4/6 installation. 

(ii) In the case of -C3 a file named mycluster.txt will 
be created containing mass, positions (a;, y, z) and veloci- 
ties (v x , v y , v z ). If you are using mcluster_sse then this file 
will also contain for each star the ZAMS masijj, the stellar 
typqj, the epoch of the staiQ its spirQ its radiuij, its lu- 
minositjlfl, its age in Myr, its metallicity (Z), its absolute V 
magnitude, its apparent V magnitudqj, B — V, its effective 
temperature (K), a random error for the V magnitude, and 
a random error for B — V . The last six are generated assum- 
ing obser vations with an 8m-class telescop e from a distance 
Rgal (see lKiipper. Mieske fc Kroupall201ll '). 

• In addition you have to specify whet her you want 
the o utput to be in iV-body units (see e.g. iHeggie fc Hutl 
120031 . option -uO) or astrophysical units (option -ul). For 
the McLuster output to serve as input for Nbody6 and 
Nbody4 this output should always be in iV-body units. 

• With the ASCII table output you can easily draw 
a colour-magnitude diagram. Use columns 18(+21) versus 
17(+20) for a diagram showing B — V versus apparent V 
magnitude (+random errors). 

• McLuster automatically computes a radial den- 
sity profile and a cumulative radial density pro- 
file. Both are by default printed to the screen. 
This may be switched off within the main rou- 
tine (parameters create_radial_prof ile = or 
create_cumulative_prof ile = 0, respectively). 



All Miscellaneous 

• The virial ratio, Q = —Eki n /E po t, where E^ in is the 
total kinetic energy of the single cluster stars and E pot their 
potential energy, can be set with the parameter -Q. Note that 
this only affects the input file for iV-body computations but 
not the stellar velocities in the table of stars (there the virial 
ratio will always be 0.5, i.e. virial equilibrium). The velocities 
get scaled within Nbody4/6 according to your choice of Q. 

• The random seed can be set with the parameter -s 
which can be any positive integer. In the case of -sO McLus- 
ter will take the local time as random seed. 



• There is an upper limit of temporary stars or orbits 
within McLuster. This number, Nmax, is set to 1.500.000 
by default, i.e. McLuster allocates memory accordingly. 
When applying mass segregation or fractality to a cluster, 
many more temporary stars/orbits have to be generated 
than finally needed. Especially if a cluster shall be mass 
segregated and fractal at the same time this number may 
easily be exceeded. In this case you should increase it in the 
main routine. 

• A few more parameters and command line arguments 
are available (see option -h and the header of the main rou- 
tine) which mostly affect flags for iV-body computations. 



A12 Examples 

If no command line arguments are passed to McLuster it 
will use the default parameter values which are specified 
and which can be changed within the main routine. By 
typing 

> mcluster 

a file test.txt is created. This default cluster has 
1000 Mg (~ 1800 stars), a half-mass radius of 0.8 pc, a 
Plummer density distribution, is in a Milky- Way tidal field 
with LSR values, uses the Kroupa IMF and has no binaries. 
The entries in the data table are in astrophysical units. With 

> mcluster -CO -uO 

the same cluster is written to the files test . input 
and test .fort . 10 but in iV-body units. This can be passed 
to Nbody4/6 as stated above. The clusters used in this 
work were created using, e.g., 

> mcluster -M100000.0 -P3 -rO . 1 -c20.0 -g2.0 -S1.0 
-CO -Gl -o R136 -fl -b0.2 -pi -s2 -Z0.01 -uO 

for the fully mass segregated JV-body model. The ar- 
guments stand for: a total mass of 100.000 Mq (-M), the 
EFF density profile (-P) with a 2d core radius of 0.1 pc 
(-r), a cut-off radius of 20 pc (-c) and a 2d power-law slope 
of -2 (-g). It is completely mass segregated (-S), the output 
is for Nbody6 (-C), and we use a GPU (-G). The output is 
named R136 (-0), we use a Kroupa IMF (-1), 20% binaries 
(-b) and ordered pairing for massive stars (-p). The random 
seed of our model is 2 (-s) and the metallicity is 0.01 (-Z). 
The output is in iV-body units (-u) since we want to pass 
it to Nbody6. 

This paper has been typeset from a TgX/ LTgX file prepared 
by the author. 



3 from SSE, in Mq 

4 from SSE, see lHurlev, Pols fc Tout]|200Cl 

5 from SSE, in Myr 

6 from SSE, in km/s 

7 from SSE, in solar units 

8 from SSE, in solar units 

9 assuming a distance Rgal from the observer which can be 
changed in the main routine 



